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Abstract. We have developed a two-dimensional orbit averaged Fok- 
ker-Planck model of stellar clusters which expands on spherically symmet- 
ric one-dimensional models to include rotation and ellipticity. Physical 
effects such as collisions, finite stellar lifetimes and bar formation [i.e., 
a non-axisymmetric component of the potential) can also be included. 
The first use of the model is to study the evolution of dense clusters 
(p ~ lO 7 M /pc 3 ) that may be expected to have existed at the centres of 
newly-forming galaxies, with the goal of verifying that angular momen- 
tum can be removed from the core of the cluster quickly enough so that 
rotation no longer prevents the formation of a massive (~ 10 2 M Q ) object. 
This could act as the seed black hole for the formation of an AGN. 



1. Introduction 

Quinlan & Shapiro (1989) developed a Fokker-Planck model of a spherically 
symmetric, dense cluster of compact stars, and found "rapid buildup of massive 
black holes in the cluster core resulting from successive binary mergers and mass 
segregation." Subsequently, they studied clusters of solar-mass stars and found 
"it is remarkably easy for massive stars to form through multiple stellar mergers 
in dense galactic nuclei" (Quinlan & Shapiro, 1990; hereafter QS). 

The goal of this project is to generalise the QS approach to two dimensions 
in order to answer the question: What about rotation? The hoped-for answer is 
that mergers and mass segregation can still produce a massive object (perhaps 
~ 1O 2 ~ 3 M ) in the core of the cluster, which could then undergo growth via 
accretion to reach supermassive size (~ 10 6_8 Mq) within a Hubble time (e.g., 
David et al, 1987) and be seen as an AGN. 



2. The Orbit-averaged Fokker-Planck equation 

Because of the non-spherical nature of our model, it is not possible to use (E, J z ) 
as our integrals of the motion as previous two-dimensional work has done (Colm 
1979). Instead we use radial (I\) and tangential (I2) action variables, which 
are adiabatic invariants as the potential evolves. Orbit-averaging is then truly 
a simple average of quantities over the 2ir change in (an) angle variable. The 
orbit-averaged Fokker-Planck equation then takes a particularly simple form: 

Wt fn = ~ W [fn (A/,)] + 28101^" (AW ^' )] - L n + G n -B n + R n 
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where f n (h,l2,t) are the distribution functions for stars of mass-type n and 
(.) are the drift and diffusion coefficients. L n ,G n are ad-hoc terms to account 
for losses and gains due to stellar mergers, while B n and R n are the analogous 
paramters for stellar evolution. Summation over repeated dummy indices i and 
j is implied. 

The coefficient calculation (partially derived by Tremaine and Weinberg, 
1984, and by Van Vleck, 1926, in different contexts) involves expanding the 
potential in action space and summing over the entire distribution / = J2n fn- 

(Aij) = -2vr 3 / dhdi 2 (V e k %f) V tj I^Wal 2 <K^i + ^2 ± e 3 n b ), 
(AiiAij) = -4vr 3 f dhdh f £ l*w 3 l 2 5{£ 1 n 1 +e 2 n 2 ±e 3 n b ). 

In the above, Q are the orbital frequencies, i.e., fl± = dE/dh, and the £ are 
integers labelling different coefficients ^e 1 £ 2 e 3 in the expansion of the potential in 
action space. Subscript b signifies "bar" , but should be understood to represent 
any individual element of the potential, e.g., that of one particular "field" star. 

The Fokker-Planck approximation remains valid as long as number N n = 
J dhdhfn 3> 1, and requires that the time step At satisfy td yn <C At <C t r . We 
have dynamical time td yn <— 10 4 y, and relaxation time t r > 10 6 y or more. 

3. The Potential 

We must iteratively solve for the self-consistent potential <E> after each Fokker- 
Planck time step. To allow for non-sphericity, we first calculate the ellipticity 
e from the overall net rotational velocity using the virial theorem (Binney & 
Tremaine, 1987). To make the problem tractable, it is assumed that isodensity 
surfaces are all of this common ellipticity. The density is thus expressed in terms 
of the homeoidal radial coordinate m? = R 2 + 

in which M n denotes the mass of star type n, and from which the new potential 
run can be calculated. This procedure is iterated until convergence is achieved, 
typically to ~ 1%, and accelerated by the Aitken "A 2 " process (Henrici, 1964). 
A final check of conservation of overall energy is also made. The m 2 grid is 
variable with each time step, with regions of larger d^/dm? and dp/ dm 2 being 
given a greater density of grid points. 

The assumption of homeiodal isodensity surfaces, along with f n , determines 
the potential $. The required integral, however, is too involved to be performed 
each time the knowledge of potential is needed in the Fokker-Planck coefficient 
calculation. After testing interpolation schemes and analytic approximations, 
only the Clutton-Brock self-consistent field method (Hernquist and Ostriker, 
1992) proved adequate in both accuracy and speed. The field method is tested 
at each timestep, and if sufficient accuracy cannot be achieved with a reasonable 
number of expansion terms, the code falls back on direct integration. 
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Figure 1. Gravitational potential $ in units of km/s as a function 
of equatorial distance r [pc] for an N = 10 6 Plummer sphere and a 
modified cluster of 1 M© stars with core radius 0.3 pc. It is probably 
wise not to put too much faith in the t = 17.2 line, as by that time 
the energy error AE/E had reached 100%. At t = 9.9, AE/E ~ 10%. 
The time unit is 0.98 x 10 6 yrs. t r (r = 0) ~ 9 x 10 6 y. 



4. The Bar 

Although some angular momentum is expected to be transported outwards via 
the effects of shear between the higher-^ inner regions and the lower-f2 outer 
regions, we are appealing to a bar-like perturbation in the potential to do most 
of the angular momentum transfer. So, should the cluster be found to have 
become unstable to formation of a stellar bar (i.e., be found to satisfy a T rot /\W\- 
style instability criterion such as that of Christodoulou et al, 1995), a non- 
axisymmetric component is incorporated into the potential. Combes & Sanders 
(1981) found that bars form over 1-2 rotational periods, and last for many more 
(> 10). This has led us to build the nonaxisymmetric potential from a fraction 
(~ 10%) of the 1 M & population that is forced to orbit in "locked step" , sharing 
a common orbital frequency and phase. We assume there is no transfer of stars 
from bar to field or vice versa. Combes & Elmegreen (1993) show that the bar 
frequency is a compromise of the orbital frequencies of its component stars, so 
our fib is set by conserving either the total energy or total angular momentum 
of those stars. Allowing the bar distribution to evolve like the field stars avoids 
any problem of inserting the bar binding energy by hand. 



5. Results 

It is standard in this game to start with a Plummer sphere distribution. To in- 
troduce rotation, we start with the rotation parameter A = J ro t \ W grav | 2 /GM t 2 ot 5 , 
and "shift" the distribution of f(J) so that total J matches desired J ro t while 
conserving number, so that the new f\(J) = for J < some J m in, and fx = 
const x / for J > J m in- A typical A value from cosmological tidal torques is 
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~0.05 (Barnes & Efstathiou, 1987). While this prescription does produce an 
overall rotation, it does so at the "expense" of the tangential velocity dispersion 
cr|, i.e., the fraction of radial orbits is enhanced and the tangential component 
of the kinetic energy is actually decreased. From a computational "proof-of- 
concept" point of view, however, this increase in radial orbits is an advantage 
as it plays into the density dependence of the gravitational relaxation. 

Figure 1 details the evolution of the potential run for N = 10 6 clusters of 
1 Mq stars. The initally Plummer-like cluster shows no change over half of a 
central relaxation time (~ 9 x 10 6 yrs). Unfortunately there was not time for a 
longer run prior to this meeting, but this demonstrates the stability of the stag- 
gered distribution/potential updating scheme. The modified cluster was run on 
a faster machine, so there was time to evolve it for ~ 2 central relaxation times. 
Note how its potential starts out deeper and deepens more rapidly than that of 
the Plummer cluster, as would be expected when radial orbits are more domi- 
nant and stars spend more time in the denser core where two-body relaxation 
has a greater effect. The corresponding densities are shown in Figure 2; the 
breaks in p(r) are caused by insufficient resolution in the innermost region.^] 



6. Conclusions and "To-do's" 

We have demonstrated the stability of this method for modeling the evolution 
of a dense star cluster using the two-dimensional Fokker-Planck equation, but 
work needs to be done to balance keeping AE/E small while maintaining a 
reasonable time step At. Ways to accomplish this include allowing a variable 
At when (AE/E) /At becomes too large (which introduces its own issues of 
energy error; these may be alleviated using the forced-time-symmetry technique 
described elsewhere in these proceedings by Piet Hut), or simply by throwing 



1 Further similar results were presented at the meeting, but in light of refinements made since, 
I will use the remainder of this space to show the equivalent newer results. 
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Figure 3. New density runs for similar clusters as in Figure 2. 

more cycles at the problem. We should also weight the small-r grid more. On 
the physical side of the problem, the initial conditions need to be chosen for 
higher v rms (a criticism of QS was that they chose initial conditions favouring 
massive object formation, but we should at least reproduce their result as a test 
case), and also not to cool of when rotation is desired. The bar can then be 
"turned on" , as can the effects of stellar mergers and evolution. 

7. Post-Denouement 

In the time between the Kingston meeting and the proceedings deadline, some 
progress towards these goals has been made. The energy error has been found 
to be largely dependent upon the differencing scheme used for the distribution 
functions f n (Ii,l2). The Chang-Cooper scheme, which in one dimension guar- 
antees that a Fokker-Planck distribution will remain positive for any At, does 
not generalise to 2-D (essentially because the calculus of extrema becomes more 
complicated). For the meeting, I used a Chang-Cooper-based method that over- 
compensated (treating all derivatives as full) but which in later tests proved 
unsatisfactory, as did undercompensating (using partials). Figure 3 shows the 
results of a simple fixed differencing weighted 90% towards the forward side, 
which allowed the 1 M case to run for more than two central relaxation times 
(i.e., over four times as long as before) before AE/E reached ~ 10%. One can 
now see the beginnings of density increase in the core. Also, the break in p(r) 
at r ~ 0.1 pc has been almost eliminated by the use of a better-optimised m 2 
grid. 

The results for clusters of the same overall mass, but with a mix of 1-, 
2- and 4-M stars is shown in Figure 4. In this case, AE/E reached ~ 20% 
after just over one core relaxation time of ~ 0.75 had passed. Here we see that 
both clusters have larger rates of central density increase than the equivalent 
all-1 M Q clusters, and that the relative increase of the rate for the cluster with 
enhanced radial orbits (quixotically labelled "rot" as described in §5) relative to 
the Plummer-like one ("norot") is much quicker than in the 1 Mq case. Both of 
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Figure 4. p(r) for clusters of 60% 1M , 30% 2M , 10% 4M by mass. 

these differences are as expected when the evolution is dominated by two-body 
gravitational relaxation. 
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